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Abstract 

We address persistence and stability of three-dimensional vortex configurations in the discrete 
nonlinear Schrodinger (NLS) equation and develop a symbolic package based on Wolfram's MATH- 
EMATICA for computations of the Lyapunov-Schmidt reduction method. The Lyapunov-Schmidt 
reduction method is a theoretical tool which enables us to study continuations and terminations of 
the discrete vortices for small coupling between lattice nodes as well as the spectral stability of the 
persistent configurations. The method was developed earlier in the context of the two-dimensional 
NLS lattice and applied to the on-site and off-site configurations (called the vortex cross and the 
vortex cell) by using semi-analytical computations [181 119) . The present treatment develops a full 
symbolic computational package which takes a desired waveform at the anti-continuum limit of 
uncoupled sites, performs a required number of Lyapunov-Schmidt reductions and outputs the pre- 
dictions on whether the configuration persists, for finite coupling, in the three-dimensional lattice 
and whether it is stable or unstable. It also provides approximations for the eigenvalues of the 
linearized stability problem. We report a number of applications of the algorithm to important 
multi-site configurations, such as the simple cube, the double cross, and the diamond. For each 
three-dimensional configuration, we identify exactly one solution, which is stable for small coupling 
between lattice nodes. 



1 Introduction 

Over the last decade, the topic of nonlinear dynamical lattices and of the coherent structures that arise 
in them has attracted considerable attention. This can be mainly attributed to a diverse set of research 
themes where the corresponding mathematical models have emerged as an appropriate description of 
the physical problem. Such areas include, but are not limited to, arrays of nonlinear-optical waveguides 
and photorefractive crystal lattices 2 , Bose-Einstein condensates (BECs) trapped in optical lattices 
(OLs) P], Josephson-j unction ladders [H, micro-mechanical models of cantilever arrays or even 
simple models of the complex dynamics of the DNA double strand £§j . These, in turn, have spurted an 
increasing mathematical interest in Hamiltonian discrete systems and produced a considerable volume 
of theoretical results, which has partially been summarized in the reviews 0. 

Arguably, the most prototypical among the discrete nonlinear Hamiltonian models is the discrete 
nonlinear Schrodinger equation (DNLS) [S], which results from a centered-difference discretization of 
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its famous continuum analog, the nonlinear Schrodinger equation (NLS) 9 j. The DNLS arises most 
notably in the nonlinear optics of waveguide arrays, where it was developed [10] as an envelope wave 
description of the electric field in each of the waveguides. It was also systematically derived as a 
tight-binding model for Bose-Einstein condensates trapped in, so-called, optical lattices In other 
settings, such as e.g. photorefractive crystals, the validity of the DNLS model is more limited, yet still a 
lot of its qualitative features can be observed even experimentally. It is for this reason that this simple 
yet rich model has been used to justify numerous experimental observations in the above-mentioned 
areas, including the formation of discrete solitons ^21; the existence of Peierls-Nabarro barriers for 
such waves |13j . the modulational instability of uniform states both in optical Jl] and in BEC jlBj 
experiments, and the formation of discrete vortices ^H] i n two-dimensions, among many others. In the 
earlier work of ^Jj and |18| I19| , we studied systematically the solutions that can be obtained in the 
one- and two-dimensional installments of the model respectively, using a combination of methods of 
Lyapunov-Schmidt reductions and perturbation theory in the, so-called, anti-continuum limit of lattice 
sites uncoupled with each other. 

In the present paper, we provide a mathematical framework and a symbolic package for the systematic 
computations of the existence and the stability of localized states in the fully three-dimensional DNLS 
model. To the best of our knowledge, this is the first time that systematic stability results are developed 
for the three-dimensional system. Earlier numerical work in j2L)l I21j revealed a considerable wealth of 
possible states, including octupoles, diamonds and vortex cubes, among many others. It is these states 
and different variants thereof that our systematic methodology addresses in the present publication, by 
formulating a symbolic approach that can be straightforwardly used to study the potential persistence 
and linear stability of any desired configuration. 

Our results on the three-dimensional DNLS model cannot be applied to dynamics of waveguide 
arrays since the third spatial direction in the optical problem represents the evolution time variable. 
However, the model is still relevant to BECs and an additional possibility for its physical realization is 
offered by a three-dimensional crystal built of microresonators |22j . 

Our article is structured as follows. Section 2 formulates the mathematical problem of interest. 
Section 3 presents the main results for the development of the computational algorithm based on the 
Lyapunov-Schmidt reduction method. Section 4 reports the application of the computational algorithm 
to three kinds of three-dimensional vortex configurations. Section 5 compares these results with direct 
numerical bifurcation results of the three-dimensional DNLS model. Section 6 concludes the paper. 
Appendix A contains typical outputs of the MATHEMATICA software package. 

2 Formulation of the problem 

We consider the discrete nonlinear Schrodinger (NLS) equation in three spatial dimensions: 

iu n + eAu n + \u n \ 2 u n = 0, n G Z 3 , t E R+, u n E C, (2-1) 

where e > is the coupling strength (which can also be thought of as the reciprocal squared lattice 
spacing) and Au n is the discrete three-dimensional Laplacian 

Au n — U n -\- ei -\- U n — ei -\- U n -\-Q2 ~\~ U<n—e2 ~t~ ^n+e3 ^n—ez ^W^, Tl G Z , 
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with {ei, 62,63} being standard unit vectors on Z 3 . The discrete NLS equation 1)2.1(1 is a Hamiltonian 
system with the Hamiltonian function 



H 



1, 



Mr, 



U 2 (Z3) 



U 4 (Z3)' 



(2.2) 



which is referred to as the energy of the discrete system. Due to invariance of the discrete NLS equation 
()2.1|1 with respect to time translation u n {t) \— > u n (t — to), Vto £ the energy H is constant in time, 
e.g. H(t) = -ff(O). Another conserved quantity is the discrete / 2 -norm, such that Q = \\u n (t)\\ 2 2 ^3^ 
and Q(t) = Q(0), which is related to the invariance of the NLS equation 1)2.1(1 with respect to the 
gauge transformation u n (t) 1— > u n (t)e te °, \/9q G K. The natural phase space of the discrete system is 

x = z 2 (z 3 )n/ 4 (z 3 ). 

Let u be an abstract vector for the triple-indexed doubly-infinite two-component sequences of real 
and imaginary parts of {it n } nG 23, such that the 2-block at the node n G Z 3 is 



u, 



Re(u r> 
lm(u r , 



Then, the discrete NLS equation 1)2.1)1 is equivalent to the dynamical system = JVi?[u], where 
H[u] : X 1— ► M 1 is the Hamiltonian function (J2.2|) and (J, V) are the standard symplectic and gradient 
operators, respectively. Operators (J, V) are block-diagonal with the 2-block at the node n G Z 3 being 



J 1 1 



1 
-1 



He(u n ) 
im(?t„) 



We consider the time-periodic (stationary) solutions of the discrete NLS equation (|2.1I) in the form 



U n (t) = CP^ 1 -^, 



n G 



(2.3) 



where all parameters have been normalized for convenience. The sequence {</>n}nez 3 i s a solution of 
the difference equation 



(1 



E = A + 6, 



n G 



(2.4) 



which is obtained by variation of the Lyapunov functional A[u] = H[u] + Q[u], such that the Euler- 
Lagrange equation VA[u]| u=( ^ = is equivalent to the difference equation (|2.4|) . 

If the localized solution cj) of the difference equation (|2,4|) exists, the time-evolution of the discrete 
NLS lattice (|2.1j) near the time-periodic space-localized solution 1)2.3(1 is defined by the linearization 



u n (t) 



J(l-6e)t 



!>„ + a n e xt + 6 n e At 



n G 



(2.5) 



where A is the spectral parameter and the sequence {(a n , b n )} ne %3 solves the linear eigenvalue problem 
for difference operators 



- 2\(p n \ 2 ) a n - cf> 2 n b n - eEa n = i\a n , -(f%a n + (l - 2\(p n \ 2 ) b n - eE6„ = -i\b n , n G 



(2.6) 



If there exists a solution with Re(A) > 0, the stationary solution 1)2.3(1 is called spectrally unstable. 
Otherwise, it is neutrally spectrally stable. 
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The difference equation (|2.4|) with e = has a general set of localized modes on n 6 Z 3 : 

(0) / e ie ", n£S, 
^ = \o, neS\ (2 ' 7) 

where S is a bounded set of nodes on the lattice n G Z 3 , S 1 - = Z 3 \S, and {0 n }nes is a set of phase 
configurations. The set {# n } n gs is arbitrary for e = 0. It is called a vortex configuration if 9 no = for 
a single node no G S C Z 3 and 6> ni / {0, 7r} for at least one m G S, ni / no- If n G {0, 7r}, Vra G 5, it 
is dubbed a discrete soliton. We address here the main questions: For what vortex configurations, the 
localized mode \2. 7| ) can 6e continued in the difference equation \2.1$ for small e > and is stable in 
the linearized discrete NLS equation 12. 6]) ? 

Although the questions and the algorithms of their solution are rather general, we would like to con- 
sider specific vortex configurations in the three-dimensional NLS lattice. Some of these configurations 
have been addressed in the previous works |20| 121 j as they are thought to be elementary building blocks 
of three-dimensional discrete structures. The three specific configurations are formulated in terms of 
the limiting solution (|2.7|) , 



(i) A simple cube consists of two adjacent planes of aligned vortex cells. More precisely, we set 
S = So® Si, where 

S l = {(0,0,0,(1,0,0,(1,1,0,(0,1,01, i = o,i, (2.8) 

such that N = dim(S) = 8. Let j be the index enumerating nodes in the contours So and Si 
according to the order (|2.8|) . Let 0y be the corresponding phase in the set {0 n }neS> where I = 0, 1 
and j = 1, 2, 3, 4. Then, the simple cube vortex configurations are 

floj— ^"^ , Oij = 8o + so 7r{j ~ 1) , j = 1,2,3,4, (2.9) 

where 6q = {0, ?[, it, and so = {+1,-1}. Two configurations with so = —1 are redundant: 
the case 0q = ir can be obtained from the case #o = by multiplication of u n by i and so does 
the case 9q = ^ from the case 9q = ^ . In what follows, we only consider the six irreducible 
vortex configurations and show that only three configurations persist for e / and only one 
configuration with 9q = ir and so = 1 is linearly stable. 

(ii) A double cross consists of two symmetric planes of aligned vortex crosses separated by an empty 
plane. More precisely, we set S = S-i Si, where 

^ = {(-i,o,0,(o,-i,0,(i,o,0,(o,i,0}, i = -M, (2-10) 

such that iV = dim(S) = 8. By using the same convention as in (i), the double cross vortex 
configurations are expressed by 

Q-ij = 7r ^ , ~ 1 \ O hj = 9o + s 7r{j ~ 1) , j = 1,2,3,4, (2.11) 

where 9q = {0, f ,vr, and so = {+1, — 1}- Similarly to the simple cube vortex configurations, 
we will consider six irreducible vortex configurations and show that three configurations persist 
for e / and only one configuration with 9q = it and sq = 1 is linearly stable. 
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(iii) A diamond consists of a quadrupole in a central plane, surrounded by two symmetric central 
off-plane nodes. More precisely, we set S = S-i © Sq © Si, where 

S = {(-1,0,0), (0,-1,0), (1,0,0), (0,1,0)}, S±i = {(0,0,±l)}, (2.12) 

such that N = dim(S') = 6. By using the same convention as in (i), the diamond vortex 
configurations are expressed by 

0oj = ttO'-1), i = 1,2,3,4, 0±i,o = 0o, (2-13) 

where 8$ = {0, f, vr, ^}. Six configurations with 0q > 8^ are redundant as they can be obtained 
from the corresponding configurations with 8^ < 8q by reflection: 8^ 8^. Three other config- 
urations with 8q = 8q = {ir, ^} and 9q = tt, 9q = ^ can be obtained from the configurations 
8q = 8q = {0, |} and 9q = 0, 9q = f by multiplication of u n by —1. One more configuration 
with 9q = and 0q = ^ can be obtained from the configuration with 8^=0 and 8q = ~ by 
complex conjugation. In what follows, we only consider the six irreducible vortex configurations 
and show that only three configurations persist for e / and only one configurations with 9q = ^ 
and 8^ = ^ is stable. 



3 Review of the Lyapunov— Schmidt reduction algorithm 



We review the main results of |18) . where the Lyapunov-Schmidt reduction method is developed to 
answer the questions outlined in Sectional 

Let O(0) be a small neighborhood of e = on R 1 . Let N = dim(5) and T be the torus on [0, 2^]^ 
for the vector 6 of phase components {8 n } n& s- We define the nonlinear vector field F(0, e) on d) G X 
and e G R 1 , such that the 2-block at the node n G Z is written by 



Fn(«M 



(3-1) 



The Jacobian D ( uF{(j),e) of the nonlinear vector field F(0, e) at the solution d> for each e G O(0) 
coincides with the linearized energy operator Ti related to the quadratic form for the Lyapunov function 
A[u] = H[u] + Q[u] such that 

A[u] = A[0] + i(V>,W</>) + 0(|Mlkx), 
where u is expressed by (j2.5l) and the 2-block of tp G X x X is defined at the node n G Z 3 by 



6„ 



The matrix operator Honi/>6lxIis not block-diagonal due to the presence of the shift operator 
E. We can still use a formal notation TL n for the "2-block" of TL at the node n G Z 3 



1 



e ( s +ei + S-ei + S+e 2 + ^-62 + s +e 3 + 8- 



e 3 J 



1 
1 



(3.2) 
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where s e .u n = u n+e - for {ei, e2, 63} G Z 3 . This notation allows us to write the matrix-vector form TLip 
in the component form TL n xjj n at each node 116Z 3 . In particular, the linear eigenvalue problem (|2.6j) 
can be rewritten in the matrix-vector form 



aHtp = iXiJ), (3.3) 

where the 2-block of a is a diagonal matrix of (1, —1) at each node n £ Z 3 . 

Let (o) = (o) (0) be the solution l(2~7|) . By explicit computations, we have F(0 (O) ,O) = and 
Ker(W(°)) = Span({e n }„ e5 ) C X x X, where = ZfyF^W.O) is block-diagonal with the 2-block 
at the node n G Z 3 given by 

(W w )»=[l Jl.ne^, (W<°))„ = 



-1 

-2i0„ 



, n £ S. 



We note that H (0) e„ = and W( )e r , 
e n at the node k £ Z 3 are given by 



(e n )i 



-2e n , n G 5 C Z 3 , where the 2-blocks of eigenvectors e n and 



with being a standard Kroneker symbol. Let PiIxIh Ker(?^ )) be an orthogonal projection 
operator to the iV-dimensional kernel of . In the explicit form, the projection operator V is 
expressed by 



Vf=(fi,f 2 )GlxI: 



(Vf) r 



1 



(e w ,f) _ 
(^nj 6 n ) 2i 



e-*(fL 



(fa), 



(3.4) 



We note that the constraint f2 = fi is added when solutions of the nonlinear vector equation F(<j>, e) = 
are considered with <j> £ X and (f> £ X. 

Since the operator Ti.^ is a self-adjoint Fredholm operator of zero index, the decomposition X x 
X = Ker(W(°)) Range (H^) is well-defined and so is the projection operator (I — V) : X x X ^ 
Ranged )). By using the Lyapunov-Schmidt reduction algorithm, we consider the decomposition 



(f) = (0) + Y / a n e n + i P £X, 

neS 

where (cp, ip) £ Range(W( )) and a n £ R for each n £ S. We note that 

neS 



(3.5) 



(3.6) 



Since the values of 6 in cj>(°\6) have not been defined yet, we can set a n = 0, Vn G <S without loss of 
generality. The splitting equations in the Lyapunov-Schmidt reduction algorithm are 

VF(4>V> (0) +cp,e) = 0, (I - 7 ? )F(0(°) (0) + y>, e) = 0. 

We note that (Z - P)H(Z - P) : Range(ft( )) i-> Range^ )) is analytic in e G O(0) and invertible 
at e = 0, while F(<f>, e) is analytic in e G O(0). By the Implicit Function Theorem for Analytic Vector 
Fields, there exists a unique solution cp £ X analytic in e G O(0) and dependent on £ T, such 
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that (f = (p(9,e) and Hvllx = 0(e) as e —* 0. As a result, there exists the nonlinear vector field 
g : T x I 1 h WL N , such that the Lyapunov-Schmidt bifurcation equations are 

g(0, e) = 7>F(0<°>(0) + <p(0, e), e) = 0. (3.7) 

By the construction, the function g(0, e) is analytic in e E 0(0) and g(0, 0) = for any 6 £ T, such 
that the Taylor series for g(0, e) in e 6 O(0) is given by 

oo 

g(0,e)=^eV fc) (0). (3.8) 
fc=l 

By the gauge symmetry, the function g(0, e) satisfies the following relation: 

VaoSM 1 , V0eT: g(0 + a oPo , e) = g(0 + a oPo , e), (3.9) 

where p = (1, 1, 1) T £ R N . The main theorem of the Lyapunov-Schmidt reduction algorithm is 
summarized as follows: 

Theorem 1 (Persistence) The configuration 4>(°\0) in \2. 7| ) can be continued to the domain e G 
0(0) i/ and only if there exists a root 0* £ T of the vector field g(0, e) in \'J. 7| ). Moreover, if the root 
0* is analytic in e £ 0(0) and 0* = 0o + O(e), the solution cf) of the difference equation \2.1$ is analytic 
in e £ 0(0), such that 

oo 

= 0(°>(0,) + ¥>(0*, e) = 0(°)(fl o ) + £ e fc (fc) (0o), (3.10) 

fc=i 

where <p^ k \6o), k £ N are independent on e. 

Implementation of the Lyapunov-Schmidt reduction algorithm is based on the analysis of the con- 
vergent Taylor series expansions (|3.8|) and (|3.10|) . Let A4 = Dgg{0,e) : R N i— > be the Jacobian 
matrix evaluated at the vector 6 £ T. By the symmetry of the shift operators S, the matrix A4 is 
symmetric. By the gauge transformation ()3.9|) . Aipo = £ W N , such that the spectrum of M includes 
a zero eigenvalue. If the zero eigenvalue of M is simple, zeros 0* of the function g(0,e) are uniquely 
continued in e modulo the gauge transformation (|3.9|) by using the Implicit Function Theorem for 
Analytic Vector Fields. 

Algorithm 1 (Persistence) Suppose that g^ k \6) = for k = 1,2, ...,k — 1 and g( K )(0) ^ for 
k>1. Let O be the root o/g (/t) (0) and = D g g^(0 ) for k > k. 

1. If Ker(M^) = Span(po) C M N , then the configuration \2. 7| ) is uniquely continued in e £ 0(0) 
modulo the gauge transformation \3.y\) . 

2. Let Ker(A^ w ) = Span( Po , P i, ...,p d J C R N with 1 < d K < N - 1 and pW : ^ Ker(A^W) 
6e the projection matrix. Then, 

(a) If g( K+1 )(0o) ^ Range(A / l^ K ^) ; i/ie configuration \2. 7| ) does not persist for any e 7^ 0. 
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(b) If g( K+1 ) {Oq) £ Range(.M( K )), the configuration \2. 7| ) is continued to the next order. Replace 

_A/f( K ) i-> p( K )_y\4( K + 1 )p( K ) 
p(«) p(«+i) : | — ► Ker(P (K) A4 (K+1) P (K) ), 

O -> o -e(^ (K) ) _1 g (K+1) (0o), 

g(M-l) g (fc+2) 

and repeat the previous two steps. 

If the algorithm stops at the order K with n < K < oo, conclude whether the vortex configuration 
continues in e or terminates at e = 0. If no K < oo exists, the algorithm does not converge in a finite 
number of iterations. 

Few remarks regarding Algorithm ^ are m place. For subsequent iterations of the algorithm, one 
needs to extend the power series 

0* = - e (m^Y 1 g (K+1) (0o) + 0(e 2 ) (3.11) 



to higher orders of e. This expansion is evaluated from the Taylor series (|3,8j) . Additionally, the 
situation when the algorithm does not converge in a finite number of iterations may imply that the 
particular vortex configuration persists beyond all orders of e and has additional parameters besides 
the parameter qo in the gauge transformation (|3.9j) . 

If the algorithm converges in a finite number of iterations, it provides a binary answer on whether 
the configuration persists beyond e / or terminates at e = 0. Simultaneously, the method enables 
us to predict spectral stability of the persistent vortex configurations. To consider stability of vortex 
configurations, we need to consider the spectrum of operators Ti and aH in the neighborhood of the 
zero eigenvalue for e € O(0). 

Using the representation (|3.2[) and the Taylor series expansion (|3.1(Jj) . we represent operator Ti by 
its Taylor series 

oo 

H = H {0) + ^e k H (k) (3.12) 
fc=i 

and consider the truncated eigenvalue problem for the spectrum of Ti: 

n (o) + eH (i) + ... + e k-i n (k-i) + e k n (k) + o( e fc + 1 )l </, = 

where /U is eigenvalue and tp G X x X is an eigenvector. 

Let a 6 Ker(M^ ) )nKer{M^ +1) )n...nKeiiM (k ~' 1) ) C R N but a Ker(M w ) for some n<k<K. 
It is clear that a has (d^—i + 1) arbitrary parameters, where = N — 1. By using the projection 
operator V in 1)3.4(1 and the relation (|3.6|) . we obtain that 

a = V ( S «« e ™ ) > ( a « e ™ ) = jD 6>^ )(0) ( o)a, 

VneS / VnGS / 



s 



where Dq4>^\6q) is the Jacobian matrix of the infinite-dimensional vector (cj>( \6) , cj)^°\6)) with 
respect to 6 T . It is clear that 

ft(0ty(0) = 0) where ip (0) = a n.e n = Dg<f>M '(0 o )a. 

Moreover, the partial (fc — l)-th sum of the power series (|3.1Uj) gives the zero of the nonlinear vector 
field Q3.1|) up to the order 0(e k ) and has (<ifc_i + 1) arbitrary parameters when 0q is shifted in the 
direction of the vector a. At the tangent space of the nonlinear vector fields (|3.1[) in the direction of 
a, the linear inhomogeneous system 

^(o)^M + w (i)^(m-i) + _ + w M^(o) = o m = l, 2, jfc - 1 

has a particular solution in the form ?// m ) = Dacfy- 171 ' (Oq)ol for m = 1,2,..., A; — 1. By extending the 
regular perturbation series for isolated zero eigenvalues of , 

^ = /) + ^(i) + ... + e ^W + 0( f W), = W e fe + 0(e fe+1 ), (3.13) 

we obtain the linear inhomogeneous equation 

H(V fe) + ... + w (fc V (0) = WcV> (0) - 

Applying the projection operator P and recalling the definition (|3.7|) . we find that the left-hand-side 
of the linear equation reduces to the form 



V 



H^DQ^-^iOo) + ... + H (fc) J D00(°)(0 o )l " = £>0g (fe) (0 o )o! = JW^a, 



while Vij)^ = a. Therefore, ^ is an eigenvalue of the Jacobian matrix and a is the corresponding 
eigenvector. The equivalence between non-zero small eigenvalues of 7i and non-zero eigenvalues of Ai 
is summarized as follows: 



Theorem 2 (Eigenvalues of TC) Let Algorithm^ converge at the K-th order and the solution <f> in 
Theorem^ persist for e ^ 0. Then, 

A^ ( J P (m ~ 1) ^ (m) J P (m_1) )e m + 0(e m+1 ) C a(H), m = k, k + 1, K, 

where X^o(Ai) and o~(TC) denote non-zero eigenvalues of matrix Ai and spectrum of operator TC, and 
p(re-l) i s an id en tity operator. 



Similarly to the computations of small non-zero eigenvalues of TC, we consider eigenvalues of the 
spectral problem (|3.3|) truncated at the k-th order approximation: 



n (0) + eH (i) + ... + e k-i n (k-i) + e k n (k) + o^+i) 



ij) = iXaij). 



By using relations e n = —iae n and TC^e n = — 2e n , Vn G 5, we can see that the linear inhomogeneous 
equation 

H(°V (0) = 2iaD e <t> { ®(0 Q )a 
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has a solution 

^ O) = ^ane n = * (o) (0 o )a, 

where ^°\6q) is the matrix extension of (j>(°\Oo), which consists of vector columns e n , n G 5*. 
Similarly, there exists a particular solution of the inhomogeneous problem 

W (o)^M + + ... + ^M^(o) = 2iaDlcj) ( - m \e )a, m = 1, 2, fc', 

in the form <p( m ) = <&( m \6o)a, where k' = (k — l)/2 if k is odd and k' = k/2 — 1 if k is even. By- 
extending the regular perturbation series for isolated zero eigenvalue of a7i^°\ 

</, = ^(0) + eV ,(D + ... + e fc-i^(fc-D + I A (V°) + ecp^ + ... + e V" ) + + 0(e fe+1 ), (3.14) 

where ^ (m) = D e ^ m \6 )oL for m = 0, 1, k - 1, ^ = & m \6 Q )a for m = 0,1,.., A;', and A = 
e fc//2 A/j/2 + 0(e fc / 2+1 ), we obtain a linear inhomogeneous problem for ^ k ' at the order 0(e k ). When k 
is odd, the linear problem takes the form 

H (o)^(k) +H (l)^(k-l) + ... + W (^(0) = i\l /2 a^ 0) . (3.15) 
When /c is even, the linear problem takes the form 

W (oty(*) + ... + iA fc/2 (yV fe ' } + ... + W (fc ' +1 V {0) ) = ^ 2 /2 V 0) - (3.16) 

By using the projection operator V, we establish the equivalence between non-zero small eigenvalues 
of the operator aTt and non-zero eigenvalues of the reduced eigenvalue problems as follows: 

Theorem 3 (Stability) Let Algorithm^ converge at the K-th order and the solution d> in Theorem 
^persist for e / 0. Let operator TL have a small eigenvalue \i of multiplicity d, such that \i = e k fi k + 
0(e fc+1 ). Then, the eigenvalue problem \'J. A]) admits (2d) small eigenvalues A, such that A = e fe ^ 2 A fc / 2 + 
0(e fc / 2+1 ), where non-zero values X/,/2 ar & found from the quadratic eigenvalue problems 

odd k: M {k) a. = ^X 2 k/2 a, (3.17) 
even k: M {k) a + ^\ k/2 C {k) a = -X 2 k/2 a, (3.18) 

where = V In^^ {6 ) + ... + n^ k ' +1 ^ W (e ) 

We note that matrix £^ must be skew-symmetric so that all eigenvalues of the quadratic eigenvalue 
problem (|3.18|) occur in pairs X k / 2 an d — X k / 2 , which is a standard feature of linearized Hamiltonian 
dynamical systems. Computations of these eigenvalues can be achieved with a simple algorithm. 

Algorithm 2 (Stability) Suppose that the solution cf) persists in Algorithm^! and compute = 
D e gW(0 ) for K <k< K . 
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1. For each order k, where eigenvalues of MS ' are non-zero, compute matrices 

2. Find roots of the determinant equation for the quadratic eigenvalue problems 

Example 1 Consider vortex configurations ([2.9)1 on the simple cube ([2.8)1 . By explicit computations, 
the bifurcation equations are non-empty at the order k = 1: 

gj$ = sm{9 u+1 - 9 U ) + sin(^ )i _ 1 - 9 ltj ) + sm{9 l+hj - 9 ltj ), 1 = 0,1, 3 = 1, 2, 3, 4, (3.19) 

where 62 j = Qoj- Roots of g*- 1 )^) occur for vortex configurations ([2.9)1 with #0 = {Oj 71 "} an d sq = 
{+1,-1}. The other vortex configurations ()2.9)1 with 9q = {^,^-} and sq = {+1,-1} terminate at 
k = 1 of Algorithm ^ The Jacobian matrix has eigenvalue \^o(M^) = 2 of multiplicity 4 for 

6q = and sq = 1, eigenvalue A^o(A4W) = —2 of multiplicity 4 for 9q = 11 and so = 1, and two 
eigenvalues X-to(Ai^) = {—2,2} of multiplicity 2 for 9 = and sq = —1. (The same result holds 
for 9q = tt and sq = — 1 by symmetry.) Algorithm ^ does not converge at k = 1, since Ker(A / f^ 1 ^) is 
four-dimensional in all cases. Our results described in Example |2] indicate that Algorithm ^ converges 
at the order K = 6 and the three vortex configurations above are uniquely continued in e G O(0) 
modulo the gauge transformation ()3.9)1 . 

4 Computations of vortex configurations 

We have created the symbolic package that performs all steps of Algorithms ^ and |2] described in 
Section |3 The user is supposed to input the configuration of active nodes S and the corresponding 
set of angles {9 n } n ^s- The package performs computations order by order to detect if the phase 
configuration persists in the Lyapunov-Schmidt reduction algorithm and if it is spectrally stable. For 
all configurations (i)-(iii) described in Section [2J we have found that Algorithm ^ terminates at the 
finite order k = K < 00 and the persistent configurations have 6* = 6q up to the order k = K, where 
6q is the root of the first non-empty correction g^(6) with k > 1. 

The package consists of two parts. The first part performs computations of the correction terms 
g^(^o) °f the vector field g(#,e) and the Jacobian matrices MS® = DQg( k \Oo) for a given vector 
6q with k < k < K according to Algorithm ^ We have confirmed that there exists K < 00 such 
that dimKer(A4^^) = 1. We have also checked that g^-^o) = 0, n < k < K for all persistent 
configurations, such that no extension of the series ([3.11)1 is necessary. All non-persistent configurations 
terminate because g' K )(#o) 7^ 0, i.e. the given vortex configuration 0q is not a root of the first non- 
empty correction g( K \0). 

Non-zero eigenvalues of M.( k ' are recorded in the first part of the package for k < k < K . By Theorem 
121 these non-zero eigenvalues give approximations of the small eigenvalues of the linearized Hamiltonian 
Ji at the order of 0(e fc ). The second part of the package performs computations of eigenvalues of the 
quadratic eigenvalue problems ()3.17)) ~ (|3.18)) according to Algorithm Finding roots of the relevant 
determinant equations is not computationally difficult because the quadratic problem ([3.17)1 is diagonal 
in A 2 and the quadratic problem ([3.18)1 is given typically by a matrix of small size. 

Example 2 Continuing Example^ we have found four non-zero eigenvalues of A^ 1 ) at each persistent 
configuration, which result in four pairs of small eigenvalues A1/2 of the quadratic problem ()3.17[) . For 
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instance, the vortex configuration (|2.9|) with 9q = ir and sq = 1 has X^o(J^i^) = —2 of multiplicity 

4 and a pair of eigenvalues A ly / 2 = ±^/2A^o(-MW) = ±2i of multiplicity 4. Continuing Algorithms ^ 

and |2 to the order k = 2, we find two non-zero eigenvalues of pWji4^P^\ which result in two pairs 
of small eigenvalues Ai of the quadratic problem ()3.18|) with a non-zero matrix £( 2 ). For the same 
selected configuration, these eigenvalues are A^o(-^2) = 2 of multiplicity 2 and a pair of eigenvalues 
Ai = ±2i of multiplicity 2. Continuing the algorithm, we have found no non-zero eigenvalues for 
matrices 

P (2) M (3,4,5) p(2) and the last 

non-zero eigenvalue for the matrix p( 2 ) P^ , such that 
K = 6. The non-zero eig envalue of p( 2 ).A/f ( 6 )p( 2 ) results in a pair of small eigenvalues A3 of the 
quadratic problem (|3.18jl . where the matrix was found to be identically zero. For the same 
selected configuration, the last non-zero eigenvalue is A^o(-M^) = —16 and a pair of simple eigenvalues 

is A3 = ± 2X^ (A4^) = ±A\/2i. The two parts of the output of the symbolic computational package 
for the selected vortex configuration are reproduced in Appendix A from the outputs of the Mathematica 
software package. 

Similar to Example |2J we have performed computations of all configurations (i)-(iii) listed in Sec- 
tion [2J Our results are summarized in Tables 1-3 where we use the following convention: the entry 
2 x 4 on the first line of Table 1 denotes the non-zero eigenvalue 2 of algebraic multiplicity 4 for a 
persistent vortex configuration. The binary conclusions on persistence and stability of the given vortex 
configuration are also listed for the reader's convenience. 



Table 1: Vortex configurations (|2.9|) on the simple cube (|2.8|) 



Si 


Persists 


Eigenvalues of H 


Eigenvalues of iaTL 


Stable 






e 


e 2 






e 


e 3 




{0, §,7T, ^} 


Yes 


{2 x 4} 


{2 x 2} 


{-16} 


{±2 x 4} 


{±2i x 2} 


{±4iV2} 


No 


{0,^,7T,f} 


Yes 


{-2 x 2,2 x 2} 


{2 x 2} 


{-16} 


{±2 x 2,±2i x 2} 


{±2 x 2} 


{±4iV2} 


No 


{f,T,f ,0} 


No 
















{^,0,4f,7r} 


No 


















Yes 


{-2 x 4} 


{2 x 2} 


{-16} 


{±2i x 4} 


{±2i x 2} 


{±4iV2} 


Yes 


{^,0,f,7T} 


No 

















Table 2: Vortex configurations (|2.11j) on the double cross (|2,10|) 



Si 


Persists 


Eigenvalues of H 


Eigenvalues of iaTL 


Stable 






e 2 


e 4 


e 






{0,f )7 r,f} 


Yes 


{-2 x 2,2 x 2} 


{-8,28 x 2} 


{±2 x 2,±2z x 2} 


{±4i,±2 v / l4x 2} 


No 


{0,f ,7T,f} 


Yes 


{-4,-2 x 3,2} 


{-8,28} 


{±2,±2i x 3,±2n/2} 


{±4i,±2\/T4} 


No 


{f,^,f ,0} 


No 












{f )0, ^,7r} 


No 












Rf,o,f} 


Yes 


{-4 x 2,-2 x 4} 


{-8} 


{±2i x 4,±2iV2 x 2} 


{±40 


Yes 


{^,0,1,71} 


No 
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Table 3: Configurations (j!H3|) on the diamond (jUTJJ) 



S-! 


Si 


Persists 


Eigenvalues of 7i 


Eigenvalues of iuTi 


Stable 








e 2 


e 4 


e 












Yes 


{-12,-6,2,2,4} 




{±2 x 2,±2 v / 2,±2i v / 3,±2iv / 6} 




No 





7T 

2 


No 















7T 


Yes 


{-2 x 2,-5± V41} 


12 


{±2i x 2, ±^-10 + 2^41, ii^lO + 2V4T} 


±2^6 


No 


7T 

2 


7T 

2 


No 












7T 

2 


TT 


No 












7T 

2 


3tt 
2 


Yes 


{-8,-2 x 3} 


-12 


{±2i x 3,±4i} 


±2i v / 6 


Yes 



5 Comparison with full eigenvalue computations 

We proceed to test the predictions of the symbolic computations against the results of direct numerical 
approximations of relevant configurations and associated eigenvalues. Starting from the anti-continuum 
limit where the solutions (|2.7|) are explicit, we use numerical continuation techniques to obtain the corre- 
sponding solutions of the difference equations (|2.4|) for finite coupling strengths of e < 0.2. Numerical 
approximations of the vortex solutions are obtained with the fixed point iterations using Newton's 
method. Once the solution {4>n} n £i 3 is obtained to the desired numerical accuracy (typically 1CP 8 ) 
on a truncated numerical domain, the eigenvalue problem ()2.6f) becomes a large matrix eigenvalue 
problem, which is fully solved using standard numerical algebra tools. The relevant eigenvalues with 
small real and imaginary parts are isolated and their dependence on e is accordingly extracted and 
compared to the theoretical predictions of Tables 1-3. The results are shown in Figures Ill9( which are 
discussed in more detail below. In general, we denote numerically computed eigenvalues by solid (blue) 
lines, while their counterparts from symbolic computations are plotted by dashed (red) lines. Pairs 
of multiple real or imaginary eigenvalues are denoted by thick solid lines, while quartets of complex 
eigenvalues are denoted by thick dash-dotted lines. 

Figure n corresponds to the simple cube configuration with Si = {0, S,7r, ^f}. As is indicated in 
Table 1, this configuration should be unstable due to an eigenvalue A ~ 2c 1 / 2 , of multiplicity four. What 
we find however is that the pair of multiple real eigenvalue splits in two identical pairs of real eigenvalues 
and a quartet of complex eigenvalues. To make things even more interesting, all four eigenvalues in 
the right-half-plane have the same real part denoted by the very thick solid line in the left panel of 
Fig. ^ The imaginary part of the quartet of complex eigenvalues is denoted by dash-dotted line in the 
right panel of the figure. If the real part of all four eigenvalues is at the order 0(e 1//2 ), the imaginary 
part of complex eigenvalues occurs at the order 0(e) with a numerical approximation A ~ 2c 1 / 2 ± lie. 
Besides these four eigenvalues in the right-half-plane, Table 1 also reports existence of a pair of double 
imaginary eigenvalues at the order 0(e) and a pair of simple imaginary eigenvalues at the order 0(e 3 ). 
All these eigenvalues are shown on the right panel of Fig. ^ by thin lines, since the pair of double 
imaginary eigenvalues splits into two pairs of simple imaginary eigenvalues. We can see that, although 
the results of Table 1 give only the leading-order approximations of small eigenvalues of the linearized 
problem, they represent adequately the pattern of unstable and neutrally stable eigenvalues. 

Figure |U describes another simple cube vortex configuration of Table 1 with S\ = {0, ^-,7r, ^}. As 
theoretically predicted, we find this configuration to be immediately unstable, due to a double pair of 
real eigenvalues at the order 0(e 1//2 ) and another double pair of real eigenvalues of the order O(e). Both 
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pairs split for small values of e but remain simple pairs of real eigenvalues for sufficiently small values of 
e. Then, a pair of the former and one of the latter collide for e ~ 0.175, leading to a quartet of complex 
eigenvalues. Another double pair of imaginary eigenvalues exists at the order 0(e 1 / 2 ) and it splits into 
simple pairs of imaginary eigenvalues. When these eigenvalues meet the continuous spectrum located 
at ±i[l, 1 + 6e], the pairs of imaginary eigenvalues generate additional quartets of complex eigenvalues 
for e > 0.113 and e > 0.125. Finally, one more pair of imaginary eigenvalues exists at the order 0(e 3 ) 
and it remains small for < e < 0.2. 

Figure El describes the third simple cube vortex configuration of Table 1 with S\ = {ir, ^-,0, 
which is linearly stable for small e. The quadruple pair of imaginary eigenvalues at the order 0(e 1 ^ 2 ) 
splits for small e into a double pair and two simple pairs of imaginary eigenvalues. All these pairs 
generate quartets of complex eigenvalues upon collision with the continuous spectrum for e > 0.1, 
e > 0.125 and e > 0.174. Therefore, the vortex configuration becomes unstable for sufficiently large 
e. A double pair of imaginary eigenvalues at the order O(e) splits for small e into simple pairs of 
imaginary eigenvalues. The additional pair of imaginary eigenvalues at the order 0(e 3 ) remains small 
for < e < 0.2. 

The double cross vortex configuration of Table 2 with Si = {0, §,7r, 4^} is shown in Figured It is 
quite interesting that both double pairs of real and imaginary eigenvalues at the order O(e) split for 
small e into simple pairs of real and imaginary eigenvalues, while the double pair of real eigenvalues at 
the order 0(e 2 ) remain double for small e. 

The double cross vortex configuration of Table 2 with S\ = {0, tt, ^} is shown in Fig. All pairs 
of real and imaginary eigenvalues are simple including the triple pair of imaginary eigenvalues at the 
order O(e) which splits for small e into three simple pairs of imaginary eigenvalues. 

Finally, the double cross vortex configuration of Table 2 with 5i = {7r, 3 ^ E ,0,|}is found to be linearly 
stable for e < 0.2 and is shown in Fig. HO The double and quadruple pairs of imaginary eigenvalues at 
the order O(e) split for small e into individual simple pairs of imaginary eigenvalues. We note that the 
domain of stability of this vortex cross configuration is wider than the one for the simple cube vortex 
configuration on Fig. |3J 

Lastly, we turn to the diamond configurations of Table 3. For the case of S-\ = and S\ = 0, shown 
in Fig. [3 the configuration is unstable due to a simple and a double pair of real eigenvalues, both 
of O(e), captured very accurately by our theoretical approximation. In addition, a complex quartet 
emerges because of the collision of a pair of imaginary eigenvalues with the continuous spectrum for 
e > 0.175. 

The second diamond configuration with S-\ = and Si = it, shown in Fig. |H1 is unstable due to 
two simple pairs of real eigenvalues, one at the order O(e) and one at the order 0(e 2 ). The simple pair 
of imaginary eigenvalues becomes a quartet of complex eigenvalues upon collision with the continuous 
spectrum for e > 0.179. The double pair of imaginary eigenvalues remains double for < e < 0.2. 

Finally, the third diamond vortex configuration with S-i = ^ and Si = 4?, shown in Fig. is 
linearly stable for < e < 0.2. In this case also, the theoretical prediction accurately captures the 
simple pair of imaginary eigenvalues at the order O(e), the triple pair of imaginary eigenvalues (splitting 
into a double pair and a simple one) at the order O(e), and the simple pair of imaginary eigenvalues 
at the order 0(e 2 ). 
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6 Conclusion 



The main benefit of the present paper is that it provides a systematic computational approach, based 
on a symbolic mathematical package to unravel the existence and linear stability of any configuration 
of interest in a three-dimensional lattice that is of paramount interest in a wide range of applications. 
The package implements the conditions of the Lyapunov- Schmidt reduction method that are necessary 
and, in the case of convergence, sufficient for persistence of relevant configurations. The reductions 
start from the anti-continuum limit of the lattice and extend the solution with respect to the coupling 
parameter e. The algorithm subsequently connects the leading-order behavior of the small eigenvalues 
of the stability problem (which bifurcate from the zero eigenvalue at e = and are responsible for the 
instability of the vortex configurations) to the Jacobian matrix of the above conditions. In so doing, it 
provides a powerful predictor that we have always found to be accurate for small values of the coupling 
parameter (typically for e < 0.1). We thus believe that this work provides a valuable tool that can be 
used to examine various configurations that may be of interest to both theoretical and experimental 
studies in optical and soft-condensed-matter systems. 
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the grants DMS-0204585, DMS-CAREER, DMS-0505663 and DMS-0619492. 



A Outputs of the computational algorithm 

A.l Output of Algorithm 1 

Jacobian matrix at the 1st order: 



Ml = 



Non-zero eigenvalues of Ml: 



Projection of M2 to Null(Ml): 



" -1 











1 








" 





-1 











1 














-1 











1 














-1 











1 


1 











-1 














1 











-1 














1 











-1 














1 











-1 _ 


-2, u[2] 


= -2, 


u 


[3] = 


-2, 


u[4] 


= -2 






1 








-1 " 






M2 







1 


-1 




















-1 


1 













-1 








1 







15 



Non-zero eigenvalues of M2: 
Projection of M3 to Null(M2): 

Projection of M4 to Null(M2): 

Projection of M5 to Null(M2): 

Projection of M6 to Null(M2): 

Non-zero eigenvalues of M6: 



u[5] = 2, u[6] = 2 



M3 



MA 



M5 


















M6 



u[7] = -16 



A. 2 Output of Algorithm 2 

Projection of iaH — XI at 1st order: 





-1-IA 2 




1 







-1-IA 2 




1 






-1 



1 X 2 







-1-|A 2 




1 



1 





-1-IA 2 






1 \ 2 

2 A 






1 





-1-IA 2 








1 





-1-IA 2 



Non-zero eigenvalues of iaH — XI at 1st order: 

w[l] = -2i, w[2] = -2i, w[3] = -2i, w[4] = -2i, w[5] = 2i, w[6] = 2i, w[7] = 2i, w[8] 
Projection of iaH — XI at 2nd order: 



1-iA 2 


-A 


A 


-1 


A 


1-±A* 




-A 


-A 


-1 


i-\x- 


A 


-1 


A 


-X 


1-iA 2 



Non-zero eigenvalues of iaH — XI at 2nd order: 

w[9] = -2i, w[10] = -2i, w[ll] = 2i, w[12] = 2i 
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Projection of iaH — XI at 3rd order: 



M3 



Projection of iaH — XI at 4th order: 

M4 



Projection of iaH — XI at 5th order: 

M5 



Projection of iaH — XI at 6th order: 

M6 = 



Non-zero eigenvalues of iaH — XI at 

W [13] = -U Sqrt[2], w[U\ = Ai Sqrt[2] 
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Figure 1: The real (left) and imaginary (right) parts of small eigenvalues of the linearized problem 
1)2.6(1 associated with the simple cube vortex configuration with Si = {0, f , vr, 4^} versus e; see the first 
paragraph of Section [5] for the meaning of the different lines. 
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Figure 5: Same as Fig. ^ but for the double cross vortex configuration with Si = {0, 4r,7r, § } 
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